#ifndef GEOMETRY_H
#define GEOMETRY_H

#include <math.h>


#define int1	signed char
#define uint1	unsigned char
#define int2	signed short
#define uint2	unsigned short
#define int4	signed long
#define uint4	unsigned long


#define PI      3.14159265358979323846264338
#define ZERO    1e-10

#define Sqr(x)      ((x)*(x))
#define Abs(x)      ((x)>=0 ? (x) : -(x))
#define Min(x,y)    (((x)>(y)) ? (y) : (x))
#define Max(x,y)    (((x)<(y)) ? (y) : (x))
#define Pos(x)      ((x)*(x>=0))                        // Positive part
#define Sgn(x)      (((x)>0) ? 1 : (((x)<0) ? -1 : 0))  // Signum

#define each(k)     for (int k=0; k<3; k++)
#define sum(a,b)    ((a)[0]*(b)[0] + (a)[1]*(b)[1] + (a)[2]*(b)[2])



// Three-dimensional vector / coordinate

class Tvec
{
  public:
    double x[3];

    // CONSTRUCTORS
    Tvec() { each(i) x[i]=0; }
    Tvec(double vec[3]) { each(i) x[i] = vec[i]; }
    explicit Tvec(double x0, double x1, double x2) { x[0] = x0; x[1] = x1; x[2] = x2; }


    // METHODS
    double Length()                             // Euclidean norm of vector
    {
        return sqrt(sum(x,x) );
    }

    bool IsNull() { return ((x[0]==0) && (x[1]==0) && (x[2]==0)); }

    void Normalize()
    {
        double length = Length();
        if (length==0) length=1;
        each(i) x[i] = x[i] / length;
    }

    Tvec Norm()                                 // Returns normalized vector
    {
        Tvec c;
        double length = Length(); if (length==0) length=1;
        each(i) c.x[i] = x[i] / length;
        return c;
    }


    void SetAvarage(Tvec source, double c)      // Weighted arithmetic mean of this and source
    {
        each(i) x[i] = (x[i]*c+source.x[i])/(c+1);
    }

    void Transform(double matrix[4][4])         // Matrix product
    {
        Tvec c;
        each(i) c.x[i] = sum(matrix[i], x);
        each(i) x[i] = c.x[i];
    }



    // OPERATORS

    double operator *(Tvec a)                   // Dot product
    {
        return sum(x, a.x);
    }

    Tvec operator ^(Tvec a)                     // Vector product
    {
        Tvec c;         // ( |a2 a3|  |a3 a1|  |a1 a2| )
                         // ( |b2 b3|, |b3 b1|, |b1 b2| )
        c.x[0] = x[1] * a.x[2] - x[2] * a.x[1];
        c.x[1] = x[2] * a.x[0] - x[0] * a.x[2];
        c.x[2] = x[0] * a.x[1] - x[1] * a.x[0];
        return c;
    }

    Tvec operator +(Tvec a)                     // Vector sum
    {
        Tvec c;
        each(i) c.x[i] = a.x[i] + x[i];
        return c;
    }
    void operator +=(Tvec a)
    {
        each(i) x[i] += a.x[i];
    }

    Tvec operator -(Tvec a)                     // Vector difference
    {
        Tvec c;
        each(i) c.x[i] = x[i] - a.x[i];
        return c;
    }
    void operator -=(Tvec a)
    {
        each(i) x[i] -= a.x[i];
    }

    Tvec operator -()                           // Inverse vector
    {
        Tvec c;
        each(i) c.x[i]=-x[i];
        return c;
    }

    Tvec operator *(double f)                   // Scalar multiplication
    {
        Tvec c;
        each(i) c.x[i] = x[i]*f;
        return c;
    }

    Tvec operator /(double f)
    {
        Tvec c;
        each(i) c.x[i] = x[i]/f;
        return c;
    }
    bool operator ==(Tvec a)                    // Vector comparisons
    {
        return ((a.x[0]==x[0])&&(a.x[1]==x[1])&&(a.x[2]==x[2]));
    }
    bool operator !=(Tvec a)
    {
        return ((a.x[0]!=x[0])||(a.x[1]!=x[1])||(a.x[2]!=x[2]));
    }
    bool operator >(Tvec a)
    {
        return ((a.x[0]<x[0])&&(a.x[1]<x[1])&&(a.x[2]<x[2]));
    }
    bool operator <(Tvec a)
    {
        return ((a.x[0]>x[0])&&(a.x[1]>x[1])&&(a.x[2]>x[2]));
    }
    double operator [](int n)
    {
        if (n<0 && n>2) return 0;
        return x[n];
    }
};



// Three-dimensional plane, defined by three points
class Tplane
{
  public:
    Tvec v[3];

    explicit Tplane(Tvec c1, Tvec c2, Tvec c3)
    {
        v[0] = c1; v[1] = c2; v[2] = c3;
    }

    Tvec GetNormal()                            // Returns plane normal
    {
        Tvec a = v[2]-v[0], b = v[2]-v[1]; // convention
        return (a^b).Norm();
    }

    double GetAngle(Tvec w, bool vectors)
    // Returns angle between plane and vector w in radians.
    // Bool "vectors" specifies whether result is in interval <0, pi> (for true), or <0, pi/2> for false
    {
        Tvec u;
        u = GetNormal();
        if (vectors)
            return (double) acos(   (w*u) / (u.Length()*w.Length()) );
            return (double) acos(Abs(w*u) / (u.Length()*w.Length()) );
    }

    void Translate(Tvec t)                      // Translate plane by given vector
    {
        each(i) v[i] = v[i] + t;
    }

    void GetEquation(Tvec &w, double& d)
    // Returns analytical equation of plane ( w[0]x + w[1]y + w[2]z + d = 0 ) via w
    {
        w = GetNormal();
        d = -1 * (sum(w.x, v[0].x));
    }

    bool ContainsPoint(Tvec t, float &a, float &b)
    // Returns whether triangle (v0, v1, v2) constains point t and
    // coordinates of t in plane in basis (v2-v0, v1-v0) via a,b
    {
        Tvec v0 = v[2] - v[0];
        Tvec v1 = v[1] - v[0];
        Tvec v2 = t - v[0];

        double bot = 1 / ((v0*v0) * (v1*v1) - (v0*v1) * (v0*v1));
        a = ((v1*v1) * (v0*v2) - (v0*v1) * (v1*v2)) * bot;
        b = ((v0*v0) * (v1*v2) - (v0*v1) * (v0*v2)) * bot;

        // X = A + a * ->AC+ b* ->AB
        return ((a >= 0) && (b >= 0) && (a + b <= 1));
    }

    Tvec Projection(Tvec p)
    // Returns projection of p to plane
    {
        Tvec v; double d;
        GetEquation(v, d);
        double t = - (v*p.x + d) / (v.Length());
        Tvec op; each(i) op.x[i] = p.x[i] + v[i]*t;
        return op;
    }
};




#define NULL_COOR           Tvec(0,0,0)
#define ZERO_COOR           Tvec(ZERO, ZERO, ZERO)
#define VecToVals(c)        (c).x[0], (c).x[1], (c).x[2]



// Three-dimensional tensor of second order

#define VEC_ROWS	1
#define VEC_COLS	2

class Ttensor
{
  public:
    double t[3][3];

    Ttensor() {  }

    Ttensor(double d00, double d01, double d02, double d10, double d11, double d12, double d20, double d21, double d22)
    {
        t[0][0] = d00; t[0][1] = d01; t[0][2] = d02; t[1][0] = d10; t[1][1] = d11; t[1][2] = d12; t[2][0] = d20; t[2][1] = d21; t[2][2] = d22;
    }
    Ttensor(Tvec v1, Tvec v2, Tvec v3, uint1 opti)		// Three vectors 
    {
		each(i) 
		{
			if (opti==VEC_COLS)
			{
				t[i][0] = v1[i];
				t[i][1] = v2[i];
				t[i][2] = v3[i];
			} else
			if (opti==VEC_ROWS)
			{
				t[0][i] = v1[i];
				t[1][i] = v2[i];
				t[2][i] = v3[i];
			} else return;
		}
    }
    Tvec Row(int n)
	{
        if (n<0 || n>2) return NULL_COOR;
        Tvec v; each(i) v.x[i] = t[n][i];
        return v;
	}
	Tvec Col(int n)
	{
        if (n<0 || n>2) return NULL_COOR;
        Tvec v; each(i) v.x[i] = t[i][n];
        return v;
	}

    Tvec operator [](int n)								// Return n-th row
    {
		return Row(n);
    }

    Ttensor operator +(Ttensor u)                       // Matrix sum
    {
        Ttensor v;
        each(i) each(j) v.t[i][j] = t[i][j] + u.t[i][j];
        return v;
    }
    Ttensor operator -(Ttensor u)                       // Matrix difference
    {
        Ttensor v;
        each(i) each(j) v.t[i][j] = t[i][j] - u.t[i][j];
        return v;
    }
    Ttensor operator *(double k)                        // Scalar multiplication
    {
        Ttensor v;
        each(i) each(j) v.t[i][j] = k * t[i][j];
        return v;
    }
    Ttensor operator *(Ttensor u)                       // Matrix product
    {
        Ttensor v;
        each(i) each(j) v.t[i][j] = t[i][0]*u.t[0][j] + t[i][1]*u.t[1][j] + t[i][2]*u.t[2][j];
        return v;
    }
    void operator =(Ttensor u)
    {
        each(i) each(j) t[i][j] = u.t[i][j];
    }

};

#define ID      	Ttensor(1,0,0,0,1,0,0,0,1)	// Identity matrix
#define Delta(i,j)	ID.t[i][j]			// Kronecker delta

Tvec Transform(Tvec v, Ttensor m);


// NUMERIC METHODS

typedef double(Ffuncclbc)(double);          // type for function f: R -> R
int Roots(Ffuncclbc *f, int max_roots, double res, double prec, double x_min, double x_max, double roots[]);
// Finds all roots of given function f.
// . res - resolution of root separation (also minimal distance of two roots)
// . prec - precision of root enumeration
// . x_min, x_max - interval of finding
// . roots - output parameter, returns root positions
// Function returns number of finded roots.

double Determinant(Ttensor t);
// Returns determinant of given matrix.
#define IsRegular(t) (Abs(Determinant(t))>ZERO)

Ttensor Inverse(Ttensor t);
// Returns inverse matrix

Ttensor Transpose(Ttensor t);
// Returns transpose matrix


bool Gauss(Ttensor t, Tvec r, Tvec &sol);
bool Gauss(Ttensor t, Ttensor r, Ttensor &sol);
// Solves system of equations (t)(sol) = (r)

int Eigenvalues(Ttensor t, double vals[]);
// Finds all eigenvalues of given matrix and returns their count.

int Eigenvectors(Ttensor t, Tvec vecs[]);
// Finds all eigenvectors of given matrix and returns their count.

#endif // GEOMETRY_H
